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A Gaussian distribution of cosmic microwave background temperature fluctuations is a generic 
prediction of inflation. Upcoming high-resolution maps of the microwave background will allow 
detailed tests of Gaussianity down to small angular scales, providing a crucial test of inflation. We 
propose Minkowski functionals as a calculational tool for testing Gaussianity and characterizing 
deviations from it. We review the mathematical formalism of Minkowski functionals of random 
fields; for Gaussian fields the functionals can be calculated exactly. We then apply the results to 
pixelized maps, giving explicit expressions for calculating the functionals from maps as well as the 
Gaussian predictions, including corrections for map boundaries, pixel noise, and pixel size and shape. 
Variances of the functionals for Gaussian distributions are derived in terms of the map correlation 
function. Applications to microwave background maps are discussed. 
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I. INTRODUCTION 



The Cosmic Microwave Background (CMB) provides the earliest obtainable direct information about the Universe. 
Upcoming experiments (see, e.g., MAP 1996; Planck 1996) will map CMB temperature and polarization fluctuations 
with unprecedented sensitivity and angular resolution. Recent analyses have demonstrated that such maps contain a 
wealth of cosmological information, allowing a high-precision determination of the fundamental cosmological param- 
eters (Jungman et al. 1996a, b; Zaldarriaga et al. 1997; Bond et al. 1997; Kamionkowski and Kosowsky 1997) and a 
clean separation of inflationary and topological defect models of structure formation (Pen et al. 1997; Hu and White 
1997). These exciting prospects will keep the CMB at the forefront of cosmology for the coming decade. 

Most theoretical and experimental results to date have focussed on the angular power spectrum of the CMB. If the 
temperature and polarization fluctuations are Gaussian, then the power spectrum contains all available information 
5^ ■ about the fluctuations; even if the fluctuations depart from Gaussian, the power spectrum is still a valuable and 
fundamental characterization of the fluctuations. For experiments on large angular scales, notably the DMR maps 
from the COBE satellite (Bennett et al. 1996), homogeneity of the Universe essentially guarantees a Gaussian 
anisotropy distribution via the central limit theorem since the experiment only probes scales outside of the causal 
horizon at last scattering (for a relevant discussion, see Scherrer and Schaefer 1995). Probes of the anisotropics 
at resolutions of one degree and smaller hold the possibility of uncovering non-Gaussian signatures. Theoretically, 
Gaussian fluctuations are well motivated, as inflation generically produces Gaussian-distributed scalar and tensor 
metric perturbations (Guth and Pi 1982; Starobinskii 1982; Hawking 1982; Bardeen et al. 1983), although this 
prediction can be circumvented in some unusual models (Allen et al. 1987; Salopek et al. 1989). Measurement of 
Gaussian fluctuation statistics at all angular scales would clear a major hurdle for inflation; conversely, discovery of 
cosmological non-Gaussianity would cast strong doubt on the inflationary scenario. If the Universe is not described by 
an inflationary model, non-Gaussianity will be an important clue in unraveling the nature of the primordial density 
perturbations. 

A variety of techniques have been applied so far to test for Gaussianity in CMB data, primarily COBE maps. 
The two-point correlation function is equivalent to the fluctuation power spectrum. In principle, the higher-order 
correlation functions provide a complete test of Gaussianity; indeed, a Gaussian distribution can be defined as one 
in which the odd (3-point, 5-point, etc.) correlations are all zero, while the even correlation functions all factorize 
into products of two-point functions. Kogut et al. (1996), following earlier analysis of the 2-year maps (Smoot et 
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al. 1994) and theoretical work (Luo and Schramm 1993; Luo 1994), have calculated certain 3-point functions for the 
COBE 4-year maps, finding the data consistent with zero (Gaussian). A drawback of analyzing N-point correlation 
functions is that computation time becomes prohibitive for higher-point functions unless restricted to specific fixed 
configurations of points, which provides a weaker test of Gaussianity than the most general case. Also, higher-point 
analysis of maps with complicated boundaries presents further difficulties. 

Another, more heuristic statistic is the two-point correlation of temperature extrema (peaks and valleys), introduced 
by Bond and Efstathiou (1987), who provide analytic approximations to the correlation function for noiseless Gaussian 
random fields. The COBE 4- year maps are also consistent with Gaussianity for this statistic (Kogut et al. 1996). 
Finally, substantial attention has been focussed on the topological genus of temperature contours (Coles, 1988; Gott 
ct al. 1990; Smoot et al. 1994; Torres et al. 1995; Colley et al. 1996; Kogut et al. 1996). The genus, as will be 
discussed below, is equivalent to one of the three Minkowski functional for a two-dimensional map. 

In this paper we concentrate on the set of statistics known as Minkowski functionals (Minkowski 1903). A general 
theorem of integral geometry states that all properties of a <i-dimensional convex set (or more generally, a finite union 
of convex sets) which satisfy translational invariance and additivity (called morphological properties) are contained in 
d+1 numerical values (Hadwiger, 1956 and 1959). For a pixelized temperature map T(n) (or, e.g., a density field), we 
consider the excursion sets of the map, defined as the set of all map pixels with value of T greater than some threshold 
v (see, e.g., Weinberg et al. 1987; Melott 1990). Then the three functionals of these excursion sets completely describe 
the morphological properties of the underlying temperature map T(n). In two and three dimensions, the Minkowski 
functionals have intuitive geometric interpretations. The key point relevant for this paper is that the Minkowski 
functionals can be calculated exactly for an underlying Gaussian field (Tomita 1990). Thus they provide a natural 
set of statistical tests for Gaussianity. Recently, Buchert and collaborators have introduced Minkowski functionals 
into the cosmological literature, primarily in the context of three-dimensional point distributions applied to galaxy 
surveys (Mecke et al. 1994; Buchert 1995; Kerscher et al. 1997). 

For a two-dimensional map, the three Minkowski functionals correspond geometrically to the total fractional area 
of the excursion set, the boundary length of the excursion set per unit area, and the Euler characteristic per unit area 
(equivalent to the topological genus). The area (Coles 1988, Ryden et al. 1989) and the genus (mentioned above) 
have been previously considered in the context of microwave background temperature maps, but not in the unified 
formal context of the Minkowski functionals. The additivity property, previously not exploited, allows an evaluation 
of the standard errors in estimating the functionals from a map, rather than resorting to Monte Carlo simulations as 
in previous work (Scaramella and Vittorio 1991; Torres et al. 1995; Kogut et al. 1996). A second property, given by 
the principle kinematic formulas, provides a recipe for incorporating arbitrary map boundaries in a straightforward 
way, a valuable advantage when analyzing upcoming maps with irregular portions of the sky cut out to eliminate 
contamination from foreground emission. 

This paper applies the formalism of Minkowski functionals to pixelized two-dimensional maps as a test of Gaussian- 
ity; the results will be directly applicable to temperature maps of the microwave background. In Sec. II, we present an 
overview of Minkowski functionals of continuous fields and explicate their relevant mathematical properties. Section 
III then considers the related case of Minkowski functionals on a lattice, or equivalently of a pixelized map. We give 
explicit formulas for variances of the functionals and corrections for boundaries. We also give explicit calculational 
algorithms for different pixelization schemes. Section IV addresses the issue of Gaussianity, giving the explicit forms 
for the Minkowski functionals for underlying Gaussian statistical distributions. Pixelization introduces corrections to 
the Gaussian functionals in the smooth case, which we derive as power series in the pixel size. Finally, we conclude 
in section V with a discussion of various issues related to analyzing maps, including noise, smoothness, and numerical 
efficiency. 



We are interested in characterizing the morphological properties of microwave background maps, the properties of 
the map which are invariant under translations and rotations and which are additive. For example, the area of the 
map which is above a certain temperature threshold is a morphological characteristic. The branch of mathematics 
known as integral geometry provides a natural tool for this characterization, known as the Minkowski functionals 
(or quermassintegrals) . In this section we give a brief review of Minkowski functionals and summarize their basic 
properties. For further details and general mathematical background, see Weil (1983) or Stoyan et al. (1987). 

Consider a convex set K in H d . The parallel set of distance r to K is the set 
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where b(k,r) is the closed ball of radius r centered at the point k. A relation called the Steiner formula can be taken 
as the definition of the Minkowski functionals 

V(K r )=^hf d )w i (Ky, (2.2) 

z=0 

where V denotes the (n-dimensional) volume. For low-dimensional spaces, the Minkowski functionals can be expressed 
simply in terms of the geometrical quantities length t{K) (for d = 1), area A(K) and boundary length L(K) (for 
d = 2), and volume V(K), surface area S(K), and mean breadth b(K) (for d = 3): 

W (K)=t(K), W 1 (K) = 2 (d=l); (2.3a) 



W Q (K)=A(K), Wx(K) = ^L(K), W 2 (K) = n (d = 2); (2.3b) 

1 27T- 47T 

W (Jf) = V(A-), Wi(JQ = -S(JQ, W 2 (K) = —b(K), W 3 (K) = — (d = 3). (2.3c) 
The Steiner formula can also be generalized to the other Minkowski functionals besides volume as 

WjiKj^^-^Wi+jiKy. (2.4) 



i=0 



The motivation for using the Minkowski functionals for characterizing morphology is the following completeness 
theorem of Hadwiger (1959). For JC d the class of convex, compact sets in R d , consider a continuous map T : JC d — > R 
which satisfies the properties of motion invariance and additivity: 

T{gK) = T(K) \/K e JC d , g G G d , (2.5) 

where Gd is the group of rigid motions in d dimensions (i.e., rotations and translations); and 

T{K X U K 2 ) + T{K X n K 2 ) = T{K X ) + T{K 2 ) (2.6) 

for K\, K 2 £ JC n with K\ U K 2 £ K. d . Then T can be expressed as a linear combination of the Minkowski functionals: 

d 

T(K) = aiWi(K), a, e R. (2.7) 

i=0 

In other words, all of the morphological information about a convex body is contained in the Minkowski functionals. 

In the following sections, we will not be concerned with characterizing a single convex set but rather a finite union 
of convex sets (i.e. map pixels), or in mathematical terms, sets which are elements of the convex ring TZ n (Hadwiger 
1956 and 1959). In two dimensions, the area and boundary length functionals have obvious generalizations to the 
convex ring, while the third functional becomes the Euler-Poincare characteristic. The Hadwiger characterization 



theorem, Eq. (2.7), generalizes to the convex ring. In this paper, we will find it convenient to work with a different 



normalization of the functionals, given by 



Mi(K) = ^Wi{K), (2.8) 



where u>d = it d / 2 /T(\ + d/2) is the volume of the unit ball in d dimensions. The three functionals of an element of 
1Z 2 are then 

M =4 Mi = ^-L, M 2 = - X , (2.9) 

where x is the Euler characteristic. 

Another set of useful relations are the so-called principle kinematic formulae (Hadwiger 1957, Santalo 1976), which 
can be written concisely as (Mecke et al. 1994) 
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/ Mi{K x n gK 2 )dg = V (*) M^^M^Ki). (2.10) 

J G d j=0 W 

The integral is over the group of motions, formally with an invariant Haar measure dg. This formula describes the 
set K\ through its intersections with a test body K 2 in random orientations. 

To analyze a map in terms of Minkowski functionals, we consider the excursion sets of the map, the map subset 
which exceeds a fixed threshold value. The threshold is treated as an independent variable on which the Minkowski 
function als o f the excursion set depend. The three functionals of interest are (up to irrelevant constant factors, 



cf. Eq. (2.3b)) the area of the excursion set A, its boundary length L, and its Euler-Poincare characteristic \ (or 
equivalently, its topological genus). The following section discusses in detail how to construct a discrete version of 
these functionals, applicable to a pixelized map. 

III. MINKOWSKI FUNCTIONALS ON A LATTICE 

The method of analysis proposed in this paper consists of comparing the values of the Minkowski functionals 
Mi calculated from an experimentally obtained map with the theoretical predictions for the expectation values and 
variances of M$ on Gaussian distributions. Carrying out this program requires calculation of the Minkowski functionals 
from given maps, properly taking into account the boundary of the observed region and the effect of pixelization. 

This section presents the necessary formalism for the analysis of Minkowski functionals on a lattice. Rather than 
treat pixelized maps as approximations to the "true" continuous temperature field, we apply the formalism of the 
Minkowski functionals directly to random functions on discrete lattices (for extensive mathematical background, see 
Serra 1982). Although some elements of this consideration are present in the literature (Hamilton et al. 1986, Coles 
1988, Likos et al. 1995), we give an independent and self-contained derivation of our results for two-dimensional 
Euclidean maps. The general expressions for the expectation values and variances of the Minkowski functionals 
are calculated for a homogeneous and isotropic random function on a regular lattice, in terms of the probability 



distributions for the random function. This will be the foundation for the analysis of Sec. IV. Explicit calculational 
algorithms for the area, boundary length, and Euler characteristic are given. We also derive formulas for the boundary 
corrections to these Minkowski functionals. 

A. General formalism 

Consider a homogeneous and isotropic scalar function / given on some regular lattice in a region S of a 2-dimensional 
plane, so that each lattice element (pixel) s is assigned a number / (s). We define the excursion set S u of the function 
/ at level u as the union of all pixels s £ S for which / (s) > u. For instance, the lattice may be a regular square 
lattice with a given step A, and the pixels s would then be squares with side A; the region S u would then consist of 
all squares s for which / (s) > u. Another possibility is a hexagonal lattice made up of regular hexagons. Defining 
the indicator function N u (s) as N u (s) = 1 if / (s) > u and N u (s) = otherwise, the excursion set is symbolically 
represented as a union 

S u = |J sN u (s) , (3.1) 

sES 

where implicitly only terms with N u (s) = 1 are present in the union. Note that a temperature map may conveniently 
be converted to dimensionless units by expressing the temperature deviation from the mean in units of the root-mean- 
square temperature deviation of the map. 

To characterize a given map /, we use the values of Minkowski functionals Mi on the excursion set S u generated 



from / at a given level u. Because of additivity of Minkowski functionals, Eq. (2.6), the following decomposition 
formula follows from Eq. (|3.l[): 



M(S U ) =J2n u (s)M(s)- N u (s 1 )N u (s 2 )M(s 1 ns 2 ) 

+ N u(si)N u (s 2 )N u { S3 )M( Sl ns2ns 3 )- (3.2) 

where the sums are taken over all different pairs, triples etc. of lattice elements s. The intersections sif^n... of pixels 
are understood in the simple geometric sense, as intersections of polygons. Note that it is only necessary to sum over 
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adjacent pixels in Eq. (|3.2|), since the Minkowski functionals M vanish on empty sets. Therefore, the series in Eq. ( |3.2| ) 
is actually finite and stops when the intersection of the lattice elements is empt y an d M (si n s 2 H ... PI s n ) — 0. For 



instance, the intersection of any two or more pixels has zero area, so the series (3.2) for the area functional Mq = A 



stops after the first term. In a square lattice, at most four squares can have non-empty intersection; in a hexagonal 



lattice, at most three hexagons intersect, and the series (3.2) stops after fewer terms. For this reason, it is generally 



simpler to an aly ze the Minkowski functionals on lattices with hexagonal symmetry. 

Equation ( pT^ ) expresses the variable M (S u ) through the function N u (s) and the constants M (s), M (s% H s 2 ), 
which are trivially calculated for any given lattice geometry (these are just the Minkowski functionals applied to 
individ ual pixels and their intersections, determined by the area and side lengths of the pixels). Direct application 
of Eq. J3.2| ) can be used to calculate the Minkowski functionals M (S u ) for a given map at a given level u by simply 
evaluating N u (s) for each pixel s. Note that despite the sums over all pairs, triples and so on appearing in Eq. ( |3.2| ), 
in fact only terms with adjacent pixels give nonzero contributions, and the required computation time is linear in the 
total number of pixels. This will also be clear from the explicit computational algorithms below. 

We now consider expectation values and variances of M (S u ). Taking expectation values of both sides of Eq. (^), 
we obtain 

(M (S u )) = ( N u (*)) M (a) - ]T (N u ( Sl ) N u (s 2 )) M ( Sl n s 2 ) + . . . (3.3) 
The averages (N u (s)), (N u (si) N u (S2)) , ... in Eq. (|3.3|) mu st b e calculated for a given random function /. Ho 



mogeneity and isotropy of / leads to a simplification of Eq. (3^3); for instance (N u (s)) will not depend on s, and 
(N u (si) N u (s 2 )) will be a function of the distance between si and s 2 only. 

Similarly, the variances of the functionals M (S u ) can be found by substituting the decomposition formula ( |3.2| ) 
into the general expression 

var [M] = (M (S u ) 2 ) - (M (S u )) 2 . (3.4) 



For instance, the variance of the area A is 



var [A] =Al\Y J ( N u N u (s 2 )) - 



(3.5) 



where A\ is the area of one pixel. The variance will in general come from two contributions, the "cosmic variance" 
due to sampling of an underlying random field, and the noise variance arising from pixel noise. We shall incorporate 
the pixel noise into the distribution for the random field (see S ec. ITv|) and treat the variance as exclusively a result of 
sampling of the underlying ensemble. The general expression Q3.4| ) can in principle be used to compute the variances 
of the Minkowski functionals analytically for given field distributions. We give a calculation for the variance of the 



area for a Gaussian distributed map in Sec. IV, along with approximations for the variances of the boundary length 
and Euler characteristic. 

A homogeneous random field may be defined by joint n-point distribution densities p n (x%, s\, s n ) for its 

values Xi on some given configurations of n points si, Si. Due to homogeneity, the densities p n depend only on 
distances between the points Si, and we write them as p\ (x), pi (x%, x%\ 7*12), P3 (xi, xi, 2:3; ri2, T23, r-&\) and so on. 
The averages of products of the indicator functions N u (s), which will frequently arise in our calculations, can be 
expressed through integrals of p n . For instance, the average (N u (s)) is related to the one-point distribution density 
Pi (x) as 

(7V„ («)) = pi (x) dx = Pi (u) . (3.6) 

J u 

Here we have introduced the cumulative distribution function P\ (u) ; it does not depend on s due to homogeneity of 
the field. Similarly, the average (N u (si) N u (s 2 )) for a pair (si, s 2 ) of points is equal to the probability that / (s) > u 
at both points and depends only on the distance r between them: 

/>oc />oo 

(N u (si) N u (s 2 )) = / dx 1 / dx 2 pi (xi,x 2 ;r) = P 2 (u;r) . (3.7) 

J u J u 

In the same manner, one can express the average of a product of the indicator functions at any n points 
(N u (si) ...N u (s n )) through the corresponding n-point distribution function P n (u; \si — Sj\), 
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(N u ( Sl )...N u {s n )) = 



dx\. 



dx n p 2 (aci, ...,x n ;rij) = P„ (u;ry) 



(3.8) 



Note that derivation of Eq. (3/7) depends on the assumption that every pi xel h as the same number of neighbors 
and so disregards the boundary pixels. The same limitation applies to Eq. (3.8). The effect of the boundary will 
be considered separately below. Also notice that that the distribution functions P n depend directly on the n-point 
distribution densities p n of the field and therefore in general cannot be reduced to the moments of the field. The 
important exception is the case of a Gaussian random field which is considered in Sec. IV. 

The following subsections consider each of the three Minkowski functionals in turn: the area, the boundary length, 
and the Euler characteristic. The value of each functional will be normalized to the total area of the map A(S). 



B. The area 



The area functional Mq = A is the simplest of the three. To calculate the value of A for a given map with N pixels, 
we count all pixels s with values f (s) > u (i.e. with N u (s) = 1) and divide by the total number of pixels: 



A(S U ) 
A(S) 



(3.9) 



This agrees with Eq. (3.2) since, as noted above, only the first term survives in that series. The expectation value of 
the area functional is then easily found: 



MS) 



Pi («) • 



The variance of the area functional is given by Eq. (3.5), which can be re- written as 



A(S V 



A(S) 



N 2 



In terms of the distribution P2 (u; r) defined by Eq. (3.7), the double sum in Eq. (3.11) becomes 

E W. N » ( fi 2)) = ^ E p 2 (u; |si - 8 2 \) ■ 



N 2 Z-^i V"" V"^-'« V^ll N 2 

Sl,S2 S1.S2 



For small lattice steps, the last sum can be approximated by an integral 



—2 E p 2{u;\si-s 2 \) 



1 



N 2 



NA 2 



P2 (it; r) 2irrdr, 



(3.10) 



(3.11) 



(3.12) 



(3.13) 



where the integral replaces the summation over steps of A, and r max is the maximum distance between points of the 
observed region. We assumed in Eq. ( 3.13| ) that every pixel S\ is surrounded by a circular region of size r max and area 
7rr max = NA 2 in which the pixels S2 is chosen, thereby disr egard ing the effect of the actual boundary of the region 
S. The existence of the boundary affects the integral in Eq. ( 3.13| ) only at large r, and it will be shown momentarily 
that the variance does not depend on the behavior at the upper limit r max , as long as the correlations between points 
become negligible at such distances. 

The variance of the area functional becomes 



A(S U ) 
A(S) 



max J 



P 2 {u; r) 2irrdr - [Pi (u) 



(3.14) 



The analytic form of the distribution function P2 {u;r) is usually difficult to obtain because the integrals in Eq. (3.7) 
cannot be evaluated. However, it is possible to draw general conclusions about the variance of A from Eq. ( 3.14 ). 
Since correlations between distant points are absent for physical reasons, the probability P2 (it; r) of two points to 
simultaneously have values of / greater than it is approximately equal to [Pi (it)] 2 at large r. If we denote by ro the 



distance at which correlations between points become insignificant so that P2 (u; r) 
integration above to obtain 



[Pi (11)] , then we can split the 



G 



A(S U ) 
A{S) 



+ f j 2vrrdr |p 2 («; r) - [Pi (u)] 2 } 
1( S .. /( ./rrdr{p 2 (u;r)-[P 1 ( U )] 2 }. 



' max W 



One can obtain an upper bound on the variance (3.15) by noting that P2 (u;r) < Pi (u), which gives 



var 



A{S) 



-MS) Pl{u) ^- p ^ u V- 



(3.15) 



(3.16) 



Although this upper bound can only serve as a rough order-of-magnitude estimate of the variance, it can be used to 
visualize the overall behavior of the variance. The peak magnitude of the variance (at u = 0) is of the order Tq/A (S) 
and is small if the size of the observed region S is much larger than the distance at which correlations become 
negligible. The same general conclusion is valid for the variance of the two other Minkowski functionals. For the case 
of Gaussian random fields, we can make further progress in evaluating the variances (see Sec. IV). 



C. The boundary length 



Now we calculate the bo und ary length functional L of the set S u , normalized to the area A (S) of the observed 
region S. We shall use Eq. (3.2), in which we only need to calculate the first two sums since the intersection of more 
than two pixels has zero length. The only intersection of two pixels that has n onz ero length is the intersection of two 
immediately adjacent pixels that have a common side. So the second sum in (|3.2[), 



E 



N u ( Sl )N u (s2)L( Sl ns2), 



(3.17) 



needs to be taken only over pairs (si,s 2 ) of adjacent pixels with a common side. The number of such immediate 
neighbors of a given pixel and the lengths of the corresponding sides depends on the lattice geometry. If we assume 
that the pixels are regular polygons with n s sides of equal length L s , then L (si fl s 2 ) = 2L S (r eca ll that the boundary 
length functional of a line segment is equal to twice the length of the segment) and the sum (3.2) becomes 



L (S u ) = n s L s J2 N u (s) - 2L S 

ses 



N u ( Sl )N u (s 2 ). 



Normalizing to the total area of the region A (S) = NA\ gives 

L{S U ) n s L s ^— \ L 
1^ Nu ( s > 



A (S) NAi 



ses 



NAi 



E N u (s)N u (s'). 



(3.18) 



(3.19) 



seS s' adj. to s 



where the second summation is perform ed ov er n s immediately adjacent neighbors s' of s, and we divide by 2 because 
of overcounting of pairs. The formula ( 3.1S ) can be directly used to calculate the boundary length functional of a 
given map at a given level u. The required computation time is aga in linear in the total number of pixels. 

Now consider the expectation value of L. We shall use Eq. ( 3.19|) and the two-point distribution function P 2 (it; r) 
i ntrod uced in Eq. (3.7). If we denote the distance between the centers of any two adjacent pixels by A, then Eq. 
(3.19) will give 



(L {S u )) _ risL s 1 
A{S) ~ A 



T^EI^M)-^ E (Nu{ S )N u { S >)) 

s(ES \ s' adj. to s 



Tl s L s 

~aT 



[P 1 {u)-P 2 {u;A)]. 



(3.20) 



Although this expression seems to depend on the lattice geometry, the ratio n s L s /Ai is the same for any regular 
lattice which contains only rt s -sided regular polygons whose centers are separated by the distance A. For instance, 
the square lattice is characterized by n s — 4, L s = A and A\ — A 2 , which gives n s L s jA\ = 4/ A; the hexagonal 
lattice with n s = 6, L s = 2A/ % /3 and A\ = A 2 V3 gives the same answer. This allows us to write the result for all 
regular lattices as 



7 



(3.21) 



Since Eq. ( p. 21 ) exhibits an explicit dependence on the lattice step A, it is important to understand its scaling 
properties at small A. It is clear that, for a continuous random field /, the distributions at nearby points are highly 
correlated, and 



Hm P 2 (n;A) = Pi(u), 



(3.22) 



which means that both the numerator and the denominator in Eq. (3.21) will tend to zero at small A. However, the 
exact value of that limit depends on the behavior of the two-point distribution function P 2 (u; A) of the random field 
at small distances A. As will be shown in Sec. IV, the limit of Eq. (3.21) as A — > for a Gaussian random field is 
finite if the random field is "smooth" as defined below. 

The general expression for the variance with help of Eq. (3.19) gives 



var 



L(S U ) 



A(S) 




N u (a) N u («') 



n s L s 



\ses 



E 

sGS 



E N * ( s ) N « (*') 

s' adj. to s I 



(3.23) 



This expression contains terms such as (N u (si) N u (s 2 ) N u (s 2 )), which are summed over all pixels s\ and s 2 as well 
as over all nearest neighbors s' 2 of s 2 . Although one can express such averages through the corresponding distribution 
P3 (m; rij), the distances b etween points depend on their relative orientation in a complicated way, which makes an 
exact evaluation of Eq. (3.23) difficult. However, an approximate expression for such terms can be obtained by noting 
that for most 3-point configurations in question the distance between s\ and s 2 is much greater than the distance A 
between the nearest neighbor pair (s 2 ,s 2 ). This suggests that we should disregard the difference between \s\ — s 2 | 
and \s\ 

w- 



and in the first approximation treat all 3-point configurations as equilateral triangles with sides \s\ — s 2 | 



s 2 | and A. Denoting the 3-point distribution P3 for such triangles by P3 (u; A, \si — s 2 |), we have approximately 

(N u ( Sl ) N u (s 2 ) N u (s' 2 )) » P 3 («; A, \ Sl - s 2 |) . (3.24) 
Similarly, the terms with products of four N u (s) can be approximated by symmetric 4-point distribution functions 

(N u (a t ) N u (s 2 ) N u ( s [) N u (s' 2 )) w P 4 (u; A, | Sl - s 2 |) , (3.25) 



and then Eq. (3.23) can be rewritten as 



var 



L(S U ) 



A(S) 



n s L s 
AiN 



E I p 2 («; r 12 ) - 2P 3 (u; A, r 12 ) + P 4 (u; A, r 12 )] 



si,s 2 
2 



[Pt(u)-P 2 (u;A)} 



(3.26) 



where n 2 = |si — s 2 |. Following the procedure and notation of Eq. (3. IS), we can approximate the double sum over 
S\,2 by an integral over r*i 2 and obtain 



var 



L(S U ) 
A(S) 



1 



[P 2 (u; r) - 2P 3 (u; A, r) + P 4 (u; A, r)} 2-nrdr 



[Pt(u)-P 2 (u;A)Y 



(3.27) 



where we have also replaced n s L s j A\ by 4/A. By arguments similar to those at the end of the previous subsection, 
the variance is again inversely proportional to the total area A (S) of the observed region. The variance of L will be 
consider ed in more d etail for Gaussian fields in Sec. |y|, where we show that the precision of the approximations made 
in Eqs. (3.24)-( 3.25| ) is sufficient to evaluate the leading term of the variance. 
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D. The Euler characteristic 



We can use Eq. (3.2) to calculate the Euler characteristic, but the expression is considerably more complicated 
because the summation is performed over all pairs, triples, etc. of pixels that have at least one common point (since 
the Euler characteristic of any non-empty simply connected set is equal to 1). For the square lattice, up to four 
pixels may have a c omm on point, whereas for the hexagonal lattice, at most three pixels intersect. A straightforward 
application of Eq. ( |3.2[ ) for a hexagonal lattice gives 



X (s u ) = j2 N u 



ses 



Y N u ( Sl )N u (s 2 ) 
Sl ns 2 #0 



E 

sins 2 ns 3 /0 



N u (si) N u (s 2 ) N u (s 3 ) 



(3.28) 



where the sums are taken over the sets of pixels which have at least one common point, that is, over the sets of 2 and 
3 adjacent pixels. The corresponding expression for the square lattice is more complicated: 



ses 



Y N u ( Sl )N u (s 2 )+ Y, N u ( Sl )N u ( S2 )N u (s 3 ) 



Y N u ( Sl )N u (s 2 )N u (s 3 )N u {s 4 ). 
sins 2 ns 3 ns 4 ^0 



(3.29) 



where the sum now also includes contributions from four adjacent pixels meeting at one vertex, along with both 
immediately and diagonally adjacent pairs. The subscript 1 is used to distinguish this definition from a slightly 
different one suggested in Refs. (Adler 1981, Coles 1988). Namely, instead of using the union of all pixels s with 
N u (s) = 1, one could take the set formed by the centers of those pixels and the lines connecting adjacent pixels (see 
Figs. la,b). As far as the topology is concerned, the only difference is that the "diagonally" adjacent pixels are now 
considered disconnected. For lattice sizes much smaller than the typical detail of the map, the Euler characteristic 
calculated for this definition of S u would be equal to that for the original definition, since in that case the diagonally 
adjacent pixels do not affect the topology of the set. The Euler characteristic of the set S u is then given by 



X2 



(S u )=YNu(s) 



ses 



Y N u ( Sl )N u (s 2 )+ Y N a (s 1 )N u (s 2 )N u (s 3 )N u (s 4 ). : 
L[sins 2 ]^o sins 2 ns 3 ns 4 ^0 



(3.30) 



with the summations goi ng ov er sets of 1, 2, and 4 adjace nt p ixels, without counting the diagonally adjacent pairs. 
(The third term in Eq. ( |3.3C| ) does not follow from Eq. (|]2) but rather is added by hand to compensate for the 
spurious "holes" inside the region hatched in Fig. lb.) 

The two definitions \i an d Xi on the square lattice are equally valid approximations of the Euler characteristic of 
the continuous excursion set, as will be any linear combination of the two. By inspection of Eqs. ( |3.29 ), (3.30) one 
finds that the simple average of X i an d X2j 



x(Su) 



Xi (S u ) +X2 (S u ) 



(3.31) 



does not contain the sum over groups of four adjacent pixels and thus is the simplest to deal with. This definition 
counts one-half of all diagonally adjacent pairs of pixels . The averaging recipe has further advantages when analyzing 
the Gau ssian c ase, s ug gesti ng that the definition ( |3.31 ) is the most suitable for square lattices. 

Eqs. (|3~28| ), (|3~29|)^(p!3l|) can be use d to c ompute the Euler characteristic for a given map on a hexagonal and a 
square lattice, respectively. As in Eq. ( |3.18 ), the sums over groups of adjacent pixels can be transformed into the 
sums over neighbors of each pixel, arriving to expressions which are manifestly linear in the total number of pixels. 

We now express the expectation value of the Euler characteristic through the distributions P n - From Eq. ( 3.3C| ) for 
the square lattice we obtain 



(X2 (Su)) 



ses 



(Nu (s)) 



£[sins 2 ]/o 



(N u ( Sl )N u ( S2 )) 



Y (Nu(si)Nu(s 2 )Nu(s 3 )Nu( Si )) 
sins 2 ns 3 ns 4 ^0 



The first term of Eq. ( 3.32| ) has already been calculated: 

Y (Nu (*)) - NPl (u) . 



(3.32) 



(3.33) 



ses 



9 




(a) (b) 



FIG. 1. Two definitions of the Euler characteristic of a random field on a lattice. The excursion set S u (the shaded curved 
region in both figures) consists of all points where the field values are above the level u. (a) The diagonally adjacent pixels 
are considered to be connected, making one connected patch with two holes (shown in thick lines and hatching); the Euler 
characteristic is 1 — 2 = — 1. (b) The diagonally adjacent pixels are not considered to be connected, which in effect reduces 
the original set S u to a smaller set shown in thick lines and hatching. The Euler characteristic is equal to 2. The difference 
between the two definitions disappears with small enough lattice size. 



In the second term, we only need to count the pairs of squares separated by the distance A. For each square there 
are four p ossible pairs, and we have to divide by 2 to prevent double counting of pairs. Therefore, the second sum in 
Eq. is 



53 (N u ( Sl )N u (s 2 )) =2NP 2 (u;A) 
£[sins 2 ]^0 



Finally, the third sum in Eq. (p.32|) is 



^ (N u ( Sl ) N u (s 2 ) N u (s 3 ) N u (s 4 )) = NPi (u; A) , 
sins 2 ns 3 ns4#0 



(3.34) 



(3.35) 



where P 4 (u; A) is the probability that / (s) > u at four points si , s 2 , S3, s 4 situated at the vertices of a square with 
side A. If we know the d istributions Pj, i = 1, 2, 4, for the random field /, then the average Euler characteristic x P er 
unit area from Eq. (3.32) is 



(X2 (Su)) 

MS) 



5 [Pi (u) - 2P 2 (u; A) + P 4 (u; A)] 



(3.36) 



The other definition ( 3.2£ ) of the Euler characteristic for the square lattice leads to the more complicated expression 

(xi (Su)) 1 



MS) 



Pi («) - 2P 2 (it; A) - 2P 2 (u; Ay/^j + 4P 3 (u; A, AV^j - P 4 (u; A) 



(3.37) 



where P3 (u; A, Ay/2) is the probability for values at th e vertices of a triangle with sides A, A, and A^/2 to be greater 
than u. Finally, the averaged definition of %, Eq. ( 3.31 ), is expressed as 
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(x(Su)) 

A(S) 



1 

A2 



Pi (u) - 2P 2 («; A) - P 2 (u; AV2) + 2P 3 (u; A, Av^) 



(3.38) 



Analogous considerations for a regular hexagonal lattice yield the following expression for the average Euler char- 
acteristic per unit area: 



(X(S U )) 2 
A(S) A 2 V3 



[Pi(u)-3P 2 («;A)+2P 3 (u;A)], 



(3.39) 



where P3 (u; A) is the probability for the field values to be greater than u at three points situated at vertices of an 
equilateral triangle with side A. 



The variance of x can be calculated in essentially the same manner as the variance of L in Eq. (3.27), with a similar 
if more cumbersome expression as the result. We omit that expression here and return to calculations of the variance 



of x m Sec. IV 



E. Boundary corrections 



A typical experimentally obtained temperature map does not cover the full sky. Moreover, certain areas of the sky 
may be excluded from the map because of the presence of foreground sources such as our Galaxy or other reasons. 
The exclusions can be represented by taking the intersection of the full sky with an appropriate "window" domai n W . 
I f we de note b y S u t he excursion set of the temperature field if it were known throughout the full sky, then Eqs. (|3.9|) , 
(3.19), (3.28)— (3.30) of the previous subsections will yield the values of the Minkowski functionals on the intersection of 
the full excursion set S u with the window W normalized to the area of W, i.e. one would obtain Mj (S u D W) /A (W). 
Although these values can be regarded as samples of the "true" functionals Mi (S u ) /A(S) normalized to the full 
sky area A(S), the existence of the boundary of W introduces systematic errors which should be removed before 
comparing the experimental values of Mi (S u ) /A(S) with theoretical predictions. This is essentially what we mean 
by boundary corrections. 

In practice, we shall not need to assume that S u is defined on a full sky S, but only that S covers a sky region 
large enough to disregard any boundary effects. Since a full-sky map is defined on a sphere, in principle additional 
complications due to curvature arise in calculations of the Minkowski functionals; one indication of trouble is the fact 
that a sphere cannot be tiled with an arbitrarily fine-grained lattice of regular polygons. However, any curvature 
corrections should be negligible as long as the pixel dimensions are small compared to the curvature scale, and any 
lattice irregularities can be dealt with individually using the general expressions Eqs. (3.2) and (3.4). For simplicity, 
we limit the following considerations to a suitably large flat region of the sky covered by a regular lattice. 

Some recipes for the boundary corrections have been proposed in the literature, in particular, for calculating the 
Euler characteristic (Coles 1988). One natural procedure (Mecke et al. 1994) is based on the kinematic property, 
Eq. ( 2.1 0| ) . Let K2 — W be a fixed-shape window domain and K\ = S th e underlying sky domain. By considering 
the measured Mi(S u n W) to be estimators for (Mi(S u n W)), Eq. fl2.10| ) can be inverted to obtain the boundary 
correction formulas 



A(S U ) = A(s u nw) 
A(S) A(W) ' 

L (S u ) _ L (s u n W) 



A(S) 
x(Su) 
A(S) 



A(W) 

x(s u nw) 

A(W) 



A(S U )L(W) 
A(S) A(W)'' 
A{S U ) X (W) 



1 L {S u ) L(W) 



A(S) A(W) 2tt A(S) A(W)' 



(3.40a) 
(3.40b) 
(3.40c) 



(note that we need to substitute L (S u ) from Eq. ( [3.40b| ) to Eq. ( 3.40c )). These formulas hold for any shape of the 
window domain W . It is of course impossible to calculate the Minkowski functionals in a larger region than that in 
which the field is actually measured, but if the measured region accurately represents the properties of the field in 
the entire region, the above boundary correction formulas will be accurate. Unless the cosmological microwave signal 
exhibits strongly non-Gaussian features on the characteristic scales of the window domain (i.e. individual pixels from 
point sources or the galactic cut) this is likely to be an excellent assumption. 
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IV. MINKOWSKI FUNCTION ALS FOR GAUSSIAN FIELDS 



In this section, we derive the expectation values of the Minkowski functionals for Gaussian fields defined on a lattice. 
First, we treat the area functional and compute its expectation value and variance. We then explore the feature of 
smoothness of random fields which is important for understanding the dependence of the Minkowski functionals on 
the pixel size. We derive the expectation value of the boundary length functional on arbitrary lattices and show 
explicitly its dependence on the lattice size and geometry. The expectation value of the Euler characteristic is then 
computed for the hexagonal lattice. 



A. Description of Gaussian fields 

For a homogeneous Gaussian field / with zero mean, the two-point correlation function 

C( Su s 2 ) = (f( Sl )f(s 2 ))=C(r), r = |*i -s a |, (4.1) 

completely characterizes the field. Customarily, Gaussian random fields are specified by their power spectrum S (fc) 
which is the Fourier transform of the correlation function, expressed in two dimensions as 

C(r) = / kJ (kr)S(k)dk. (4.2) 
Jo 

However, our present considerations are based on the physical space (r) rather than on the momentum space (fc), so 
the correlation function C (r) is more relevant in this context than the power spectrum. 

For simplicity, we assume that the variance of the field (which is the same at all points) is (/ 2 ) = C (0) = 1. (This 
can always be achieved by normalizing the field, / = fa, if the original field had variance a 2 . To recover the original 
variables, one only needs to divide the function values / by a and the correlation function C by a 2 in all formulas.) 
The correlation function of the normalized field will then satisfy \C (r)\ < 1. Also, since correlations between distant 
points should be absent, the correlation function must satisfy C (r = oo) = 0. 

Calculations of the Minkowski functionals require the distribution functions Pi (it), Pi (it; r), and so on, introduced 
in the previous section. These functions describe the probabilities of encountering field values smaller than the level 
it simultaneously at some given points si, s%, s n and c an b e expressed through the n-point probability densities 
Pn (fi, ■ fn) for values f t = f (s,) at these points, cf. Eq. (ET~ 



DC 



P n (it) = / dh... / df nPn (fx, /„) . (4.3) 

J u J u 

Since the field / is Gaussian, all its n-point densities are finite-dimensional Gaussian distributions with probability 
densities of the form 

MA M = { 2«r / ^ eM - B " ,i,i) (44) 

with suitable coefficients B % i forming a positive-definite matrix. As is well known, the matrix B^ is the inverse of the 
correlation matrix Cy = (f (sj) / (sj)). Since the field is homogeneous, the coefficients dj are completely determined 
by the separations between points s, and Sj, i.e. Cy = C (\si — Sj\). Note that the presence of Gaussian pixel noise 
can be easily described at this stage. Imposition of noise with fixed standard deviation a 2 on a Gaussian field with 
given C (r) increases the value C (0) by a 2 while not changing C (r) at any other point. (For a discrete map, the 
function C (r) is only available on a discrete set of values r = 0, A, and so on.) We give the explicit form of P2 (fi, fi) 
in Appendix |A], including the pixel noise correction. 



B. Area and its variance 



We now consider the area functional of a Gaussian field. As follows from Eq. ( 3.10| ) , the expectation value of the 
normalized area functional is equal to the value of the distribution function P\ (u) defined in the preceding section. 
The one-point probability density for the values of the field is 



12 



Pi (/) 



and therefore the one-point distribution function is 



/2vr 



exp 



Z! 

2 



Pi(u) 



Pi (/) df 



erf 7i) 



(4.5) 



(4.6) 



This is the expectation value of the area of the excursion set S u per unit a rea of the region. 
The variance of the area functional is approximately given by Eq. ( |3.15 ), 



A(S U ) 
A(S) 



A(S) 



2irrdr 



(4.7) 



Since the integrand is assumed to vanish at r > r , the upper limit in Eq. (4.7) can be extended to infinity. As shown 
in Eq. (B13) of Appendix |b| the two-point distribution function Pi (u;r) can be expressed through the correlation 
function C (r) as 



P 2 (u-r) = [Pi(u)f + 



C(r) 



dC 



2ttV1 - C 2 



Then we can transform Eq. (4.7) into 



A(S U 



A(S) 



MS) Jo 27rrdr Jo 



exp 



dC 



I 



r 2 (C) dC 
A(S)J 2VT~~C 2 



2nVT~C 2 

„,2 



1 + C 



exp 



(4.8) 



1 + C 



exp 



1 + C 



(4.9) 



wher e w e have changed the order of integration and used the inverse r (C) of the correlation function C (r). Thus, 
Eq. (4.9) allows one to calculate (at least numerically) the variance of the area functional at a given level u for any 
given correlation function C (r). The maximum of the variance (at u = 0), as can be seen from Eq. (^^) and as we 
have checked for some sample functions C(r), is proportional t o r 2 / A (S), where r is the typical distance at which 
correlations vanish, in agreement with the qualitative estimate ( 3.16 ). 

W e can also obtain the asymptotic form of the variance in the limit of large levels \u\, for which the integral in Eq. 
(4.9) is dominated by the neighborhood of C = 1 (small distances). We assume a Taylor expansion of the correlation 
function C (r) at small distances, 



C (A) = 1 - Ci 



O (A p+1 ) 



(4.10) 



where Ci, r\ and p < 2 are some constants. Then by expanding the integrand in (1 — C) we obtain the leading term 

-(4+p)/2p 



A(S V 



A(S) 



A (S) C{ 



2/p 



V 



2p 



exp(- y 



(4.11) 



The square root of the variance ( 4.11 ) is proportional to exp (— u 2 /4) and at large enough u becomes large relative 



to the "signal" itself which, as follows from Eq. (4^6), is proportional to exp (— it 2 /2). This is to be expected for all 



Minkowski functional, since the signal at large u comes from very few map pixels. Thus, we can only use features of 
the Minkowski functional profiles at levels u within a few standard deviations from the average. 



C. Smooth random fields 



A random field is continuous if the values at nearby points are highly correlated. For a continuous field, the limit 
of the correlation function at vanishing distance must be 

HmC(A)=C(0) = l, (4.12) 
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i.e. the correlation function must be continuous at A = 0. We shall call a continuous random field smooth (with 
a somewhat relaxed rigor) if its spatial derivatives are finite, i.e. have finite variances. To see what this means 
mathematically, consider a homogeneous one-dimensional random field / (x) with a given correlation function C (r) . 
The variance of a finite-differenced derivative (/ (x + A) — / (x)) /A is 



A) - /(*)]' 



l-C(A) 



(4.13) 



We call the field / smooth if Eq. (4.13) has a finite limit at A 
function is regular at A = with vanishing first derivative, 



0. This happens for instance if the correlation 



C(A) = 1 



ic"(0)A 2 



C(A 3 ) . 



(4.14) 



In general, if the asymptotic form of the correlation function for small A is a power law, Eq. ( 4. 1 C| ) , then if p < 2 the 
field is not smooth. Physically, a random field is smooth if its small- wavelength modes are suppressed strongly enough 
so that the field values at nearby points are highly correlated. One can show, for instance, that if the small-wavelength 
modes in the power spectrum of a Gaussian random field (cf. Eq. (4.2)) are exponentially suppressed, as in the power 
spectrum 



S (k) cx k n exp (— ak) as k — > +c 



(4.15) 



then the correlation function C (r) satisfies Eq. ( |4.14 ) and the field is smooth. 

Although the CMB temperature is expected to be smooth because of Silk damping of CMB fluctuations on scales 
below a few arcminutes, the behavior of t he co rrelation function on larger scales, particularly the pixel size of a given 
experimental map, may differ from Eq. ( 4.14 ). This could make the CMB temperature field effectively non-smooth 
on the pixelization scale. As we show below, excursion sets of non-smooth fields have a fractal shape manifested by 
the scaling of the L and x functionals with the pixel size. In our subsequent analysis, we shall not assume that the 
correlation function possesses an expansion of type (4.14) on the relevant pixel scales. 



D. Boundary length for a Gaussian field 



The average boundary length L (S u ) per unit area can be found using Eq. ( 3.21 ). We need to know the prob ability 
distribution P2 (u; A) for v alues /i^ at two points S1.2 separated by the pixel size A. Substituting Eq. (Bll) from 
Appendix ||] into Eq. ( |3.2l| ), we obtain 



(L(S U )) 

MS) 



2 

7rA 



dC 



'c(A) v 7 ! - C 2 



exp - 



u 



l + C 



(4.16) 



For small enough A the values of the correlation function C (A) are close to 1, so we can expand the above integral 
in powers of the small parameter e = 1 — C (A) and obtain 



(L (S u )> 2V2e 



A(S) 



7rA 



exp 



15 



u 
45u 2 



1 



1 + 
15u 4 



12 



3 - 6u 2 + u 

160 
105 - 420u 2 



210u 4 - 28u 6 



2688 



55296 



O U b 



(4.17) 



This expansion holds, as does Eq. ( [3.21 ), for all regular lattices; the expansion parameter e is typically a power of the 
lattice step A. We note that the depend ence on A appears already in the leading term. If the fi eld is smooth and its 
correlation function is described by Eq. (4.14), then it is easily seen that the leading term of Eq. (4.17) is independent 
of A in the limit of small A. However, for non-smooth fields the dependence on A does not go away and may lead 
to formal divergence of L at small A. For instance, if the correlation function satisfies Eq. ( 4.10| ) in some range of 
A, then the leading term scales with A as A p / 2 ~ x . The scaling property of the leading term can be interpreted as a 
manifestation of fractal geometry of the conto ur lines with fractal dimension 2 — p/2 > 1. 

It is interesting to compare the series ( 4.17j ) and the result for the length of level contours of a continuous smooth 
Gaussian field in two dimensions [Adler 1981], 



(L(S U )) _ y/^(0) 

A (S) 2 



exp 



u 

~2 



(4.18) 
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With the correlation function (4.14), the leading term of Eq. ( 4.17| ) reproduces the continuous limit ( 4.18| ) up to the 
factor 4/7r. This discrepancy is due to the fact that the length of a discrete approximation of a line by a fixed regular 
lattice differs from the length of the line by a geometric factor. This factor is easiest to derive for the square lattice, 
where a curved line is approximated by horizontal and vertical straight line segments. For lattice steps A much 
smaller than the typical radius of curvature, the curve is locally well-approximated by a straight line. If a straight 
line segment of length L makes an angle a with the lattice (see Fig. 2a), then the length of its lattice approximation 
is L = L (cosa + sina). Averaging this over all angles a (it suffices to consider < a < 7r/2), we obtain 



tt/2 4 

L (cos a + sin a) da = —L. 



(4.19) 



This means that the lattice approximation of the length of a curve with isotropically distributed tangent angle is 
wrong by a factor of 4/7r on the average. By a similar argument it can be shown that any tiling by identical regular 
polygons gives rise to the same factor A/n. (In the ca se of non- regular tilings, the correction factor will generally 
differ from 4/7r, and also the higher terms of the series (4.17) will be different.) 

A more intuitive way to obtain this correction factor A/tt is to consider a square with an inscribed circle. The 
length of the lattice approximation to the circle is equal to the perimeter of the square (see Fig. 2b). The ratio of the 
circumferences is A/ir. The circle effectively averages over all orientations of the line, therefore the result also holds 
on the average for curved lines with isotropically distributed tangent angles. 




Lcosoc 

(a) (b) 



FIG. 2. Geometric correction factor for curves on a lattice, (a) Approximation of a straight line segment L making an angle 
a with the lattice has length L (cos a + sina). Averaging this over the angles a yields 4L/n. (b) The circumference of a circle, 
2nR, is approximated by a rectangular lattice. The result, even for an arbitrarily fine lattice, is 8-R, which is off by the same 
factor 4/-7T. 



Apart from the geometric correction factor, a comparison of the series (4.17) with the continuous limit (4.18) shows 
that higher terms o f the series contain polynomials in u which distort the Gauss ian profile. The resulting small 
deviation from Eq. ( 4. lg ) is a direct consequence of pixclization. While Eq. ( 4.17 ) holds for any regular lattice, in 
particular for square and hexagonal lattices, in asymmetric lattices the coefficients of the series for L depend on the 
pixel geometry. 



Application of the general formula of Eq. (3.27) for the variance to the Gaussian case requires calculation of the 
distribution P4 (u; A,r) for the two pairs of closest neighbor points separated by the distance r. In Appendix [b|, an 
expansion of the distribution P4 (u; A, r) in powers of e = 1 — C (A) is found. We retain only the leading term of this 
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expansion, because the error of the approximation in Eq. (3.27) is typically worse than the contribution O (e) of the 
higher terms of the expansion. The resulting expression for the variance of L is 



var 



L(S U ) 

MS) 



2 

7rA 



2c 



(l-C(r) 2 ) 



2N-1/2 



exp 



1 + C(r) 



exp l—u 



2-nrdr. 



(4.20) 



As in the previous section, one can transform this integral to the variable C (r) and show that the asymptotic form 
of the variance (4.2C) at large |u| is proportional to exp (— u 2 /2). 



E. The Euler characteristic 



Finally we consider the Euler characteristic of the excursion set per unit area. We first cons ider the simpler case of 
a hexagonal lattice for which the Euler characteristic is given by Eq. ( |3.39| ). Using Eq. ( p28| ) from Appendix || and 
expanding it in powers of the small parameter e = 1 — C (A), we obtain the series 



(X(S U )) 27iexp(- U 2 /2) 



A(S) 



(2tt) 3/2 A 2 



1 



-83- 



105 - 105w 2 



36 

21u 4 - M 6 



7- 



15 - 10u 2 



540 
945 - 1260w 2 



90720 



-319- 



378u 4 - 36u 6 



6123600 



(4.21) 



The leading term of the series (4.21) coincides with the formula of Adler (1981) for the average Euler characteristic 
of excursion sets of a smooth Gaussian field, 



(X) 
A(S) 



\C"(0)\ 



(4.22) 



The case of a s quare l attice is more complicated. We can use the methods of Appendix |b] to obtain a series 
expansion of Eqs. |3~3^)-(^38|). However, since the point configurations for the required distributions P3 and P4 
contain neighbors separated by distances A and A-s/2, there are now two small parameters, namely e = 1 — C (A) 
and be = 1 — C ( A-\/2) , where we have introduced the parameter 6 which is of order 1 and which characterizes the 
behavior of the correlation function C (r) at smallest distances available on the lattice, 



1 - C (AV2) 
" l-C(A) 



(4.23) 



By definition, b > 0; moreover, one can derive the condition b < 2 from the requirement of normalizability of the 
four-point distribution for the corners of a closest neighbo r squ are (cf. Appendix |A|, Eq. ( |A3[ ) with the correlation 
matrix of Eq. (A6)). A smooth Gaussian field satisfying Eq. ( 4.14 ) is described by b « 2, and the deviation of b from 
2 is of order 1 — C (A). From a given experimental map, the parameter b can be estimated in linear time. If the value 
of b is significantly smaller than 2, namely if 2 — 6^1 — G (A), we conclude that the underlying field is not smooth 
on scales A. 



After computing the expansions of Eqs. (3.36)— (3.3S) in e (while keeping the parameter b constant), one finds in 
each case an expression similar to Eq. ( 4.21 ), except that both the leading term and the coefficients of the polynomials 
in u in t he hi gher t erms depend on b. However, the most significant difference is that the leading term of the expansion 
of Eqs. ( 3.36 ) and ( [3.37| ) is not of order e but of or der y/ e, and also the expansion parameter of the series is y/e instead 
of e. The first two terms of the expansion of Eq. (3.37) are 



(X2 (S u )) = 1 
A(S) V2tt 2 A 2 



2 arccos ■ 



4-6 



Vb 



36 



arccos 



4-6 



exp 



u 
~2 



yb(^b) 
(27Fw uexp 



u 
~2 



(4.24) 



At 6 = 2 — O (e) , which corresponds to a smooth field, the first term as well as all terms with non- integer powers 
of e vanish, and the series (4.24) beco mes a gain a series in e. For non-smooth fields, however, Eq. (4.24) does not 
correspond to its continuous limit Eq. ( 4.22 ). 

A similar res ult is obtained for xii El- (3.36). The situation is however different with the Euler characteristic 
defined by Eq. (3.38): all terms with non-integer powers of e cancel, and the result is 
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(x(Su)) _ VH^b) ( u 

A (5) " (2^)3/2 A 2 P ^ 2 1 



(3 - u 2 ) (4 + 26 - 6 2 ) 



12 (4 - b) 

(b 4 - 6b 3 + 6b 2 + 2b+ 18) (u 4 - Wu 2 + 15) 

180 (4 -bf 
(6 6 — 10b 5 + 306 4 - 206 3 - 126 2 



- 80) (u e - 2lu 4 + 105u 2 - 105) 



3360 (4 - by 



+ 0(e 5 ) (4.25) 



The series ( 4.25| ) is reduc ed to Eq. ( 4.21 ) at b = 1. The leading term of this series is, up to a constant, the same as 
the continuous limit Eq. ( | 4 . 2 2| ) ; the correspondence is exact at b = 2. For this reason, we suggest that the "averaged" 
definition of x^ Eq. (3.38|), be used for square lattice calculations. As in the previous section, we note the scaling 



properties of Eqs. ( 4.21 )-( 4.25) with A and interpret them as a manifestation of the fractal structure of the excursion 
set. For smooth fields, the dependence on A cancels out, while for non-smooth fields it may lead to a divergence of \ 
at small A. 

Th e vari an ce of \ can be estimated similarly to the variance of the boundary length L, except that the expressions of 
Eqs. ( |3~3"8l) , (|3~39|) for x contain 3-point distributions P3 and therefore its variance requires calculation of the 6-point 
distribution Pq. This calculation can also be performed by methods presented in Appendix ^ and, as in the case of 
the variance of L, only the leading term of the resulting expansion suffices up to the precision of approximation. We 
give the result for the variance of the Euler characteristic only in the simplest case on hexagonal lattice: 



x(Su) 
A(S) 



^ " 'max JO 



u 2 {l-C) + C(l + C) 

- c) (1 + cy 



■ exp 



1 + C 



u 2 exp (— u 2 ) 



rdr, 



(4.26) 



where we imply C = C (r). The properties of the variance of x ar e similar to those of the variances of the other two 
Minkowski functionals, as has been discussed above. 



V. SUMMARY AND DISCUSSION 



In this paper, we have proposed the Minkowski functionals - area, boundary length, and Euler characteristic - of 
excursion sets as a probe of a map's Gaussianity. They are all linear functions of the number of map pixels, and thus 
are easy to compute from a map. Gaussian pixel noise is straightforward to include, and irregular map boundaries, 
arising from partial sky coverage or cuts of the data, can be accounted for in a natural manner. Finally, for Gaussian 
distributions the functionals can be calculated exactly and their variances estimated given the correlation function of 
the distribution, making a test of Gaussianity straightforward and eliminating the need for Monte Carlo simulations. 

We have presented a complete formulation of Minkowski functionals on a two-dimensional pixclizcd map, including 
explicit boundary corrections and pixel-size dependences. We propose an alternative definition of the Euler charac- 
teristic for pixelized maps which possesses nice calculational properties compared with previously used definitions. 
In addition, we provide explicit forms of the functionals and their variances for the case of Gaussian distributions, 
including pixelization corrections. 

Minkowski functionals are simplest to calculate for maps with hexagonal pixelation schemes, because no more 
than three pixels can ever adjoin a single point and all adjoining pixels are the same distance apart. This property 
also makes estimation of the Euler functional more straightforward and allows easier estimation of the variances 
of Gaussian map functionals. Hexagonal pixelizations also have other nice properties like optimal smoothness and 
regularity (Tegmark 1996). Analysis of maps with square pixels will be a bit more technically involved, but not 
fundamentally different. Actual full-sky maps will always have irregular points in the pixelization lattice (since the 
sphere cannot be tiled by identical polygons) which will need to be handled on an individual basis. No conceptual 
difficulties arise in extending the Minkowski functionals to a sphere as long as the pixelization scale is small compared 
to the curvature scale, which will always be the case for maps of moderately high resolution. 

It is impossible to make a general statement about how well these functionals will distinguish non-Gaussianity 
in a map because the variety of non-Gaussianity is endless. Since Minkowski functionals are ensemble or map 
averages, it is likely they will be best at picking out non-Gaussian distributions which are spread over an entire 
map, as opposed to isolated, sharp features which will be lost in the averaging. Even if the primordial perturbations 
are Gaussian-distributed (as expected in inflationary models), potentially interesting non-Gaussianity might arise 
from weak gravitational lensing (Bernardeau 1996) or from Sunyaev-Zeldovich distortions. For a specific kind of 
non-Gaussian contribution, simulations of the resulting sky will be required to determine the utility of Minkowski 
functional analysis. 
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Any test of primordial non-Gaussianity in CMB maps is likely to be challenging with real data. Just as noise 
correlations can mimic CMB power (Dodelson and Kosowsky, 1995), any noise correlation will create non-Gaussianity 
in a map (Kogut et al. 1994). A meaningful Gaussianity test requires detailed understanding of an experiment's 
noise properties. The best bet for Gaussianity tests are those experiments with simple noise properties; the MAP 
satellite design, for example, was driven largely by the desire for very simple noise properties. Other unavoidable noise 
correlations are induced by projecting out foreground contamination, but these correlations can likely be understood in 
detail through modelling and simulations. Conversely, for experiments with complicated noise properties, Minkowski 
functional analysis may aid in characterizing the noise as long as the underlying fluctuations are close to Gaussian. 

Gaussian-distributed initial fluctuations are a decisive prediction of inflationary models. Upcoming CMB data sets 
in the form of high-resolution maps will allow meaningful tests of this prediction. The theoretical framework presented 
in this paper refines and extends the current menu of Gaussianity tests, and we anticipate that Minkowski functionals 
will be a basic component of future microwave background map analysis. 
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APPENDIX A: GAUSSIAN PROBABILITY DENSITIES IN SEVERAL VARIABLES 



We first consider the simplest case of the Gaussian two-point probability density pi (x, y; A) for field values x, y at 
points si,2 separated by the distance A. For a normalized Gaussian field, the correlations are 

<* 2 > = (v 2 ) = i, 

(xy) =C(A). 

We can easily diagonalize the correlation matrix if we notice that 

((x + y){x-y)) =0, 

'(x±yf) =2(1±C(A)), 



which means that (x + y) and (x — y) are independent Gaussian variables with known variances. One can then write 
the joint distribution for (x ± y) as a product of two Gaussians, 



P2 (x, y; A) dxdy 



dxdy 



: exp 



2Wl-C(Af 



1 / (x + yf (x-yf 



4 I 1 + C (A) l-C(A) 



(Al) 



More generally, a multivariate Gaussian distribution in n variables x\, x n is completely specified by a symmetric 
correlation matrix = (xiXj), and the joint probability density is expressed through the inverse matrix C _1 by 



p n (x;Cij)d n x 



: exp 



1 - 



'■J=l 



(A2) 



If the matrix C,j is explicitly diagonalized and its eigenvalues A^ and the corresponding orthonormal eigenve ctor s E J , 
are known, then the probability density p n (x; CV,) can be written as a product of Gaussians similar to Eq. ( Al), 



Pn ( x ! Cij) 



(2ir) n/2 v /detQ 



: exp 



1 ™ 



FJ T ■ 



(A3) 



However, this explicit form of the distribution density will not be as useful as its representation as a Fourier transform 
of its generating function p n (k; Cy), 
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p n (x;Cij) = J j^ytPnfcCij) = J ^yr ex P 



1 - 

«kx— — djkikj 



i,3=l 



(A4) 



This for m o f the probability density will be used in Appendix |b] for calculations of the probability integrals P n defined 
by Eq. ©. 

The distributions for values of a homogeneous and isotropic Gaussian field at the vertices of an equilateral triangle or 
a square are found as special cases of the multivariate Gaussian distribution. The three-point distribution p% (x, y, z; C) 
is characterized by a single parameter, the correlation C = C\% — (xy) = (yz) — (zx) between any two different points, 
while Cu = 1. The correlation matrix Cij is 




(A5) 



Similarly, one obtains the four-point distribution p^ (w, x, y, z) for values at the vertices of a square. This time the 
distribution has two free parameters, namely the correlations between the adjacent points C\i = C and between the 
diagonally opposing points C23 = D. The correlation matrix Cy in that case is 



( 1 C C D 

C 1 D C 

C D 1 C 

\D C C 1 



(A6) 



The explicit forms of the distributions P3 (x, y, z; C) and p^ [w, x, y, z\ C, D) can be obtained using Eq. (A3), but we 
omit them as they will not be useful for our calculations. 

Finally we briefly describe how a known level of Gaussian pixel noise present in experimental maps can be represented 
in this formalism. If we assume that Eq. (A3) describes the n-point distribution obtained from a noiseless map, then 
pixel noise with a fixed variance a 2 (in units of the field variance) will simply increase the diagonal elements of the 
matrix Cij by <r 2 , while not changing any other correlations. Thus pixel noise can be straightforwardly incorporated 
by making the replacements Cij —> Cij /(l + er 2 ) for i ^ j in the correlation matrix. 



APPENDIX B: THE GAUSSIAN PROBABILITY INTEGRALS 



In this appendix we give some derivations for the distributions P n (u) which define the probability for the values of 
a Gauss ian random field at n given points to be above the threshol d va lue u. Generally, these distributions, as defined 
by Eq. (4.5), involve multiple integrals of Gaussians of the type (A2) and cannot be exactly evaluated. However, if 
the correlation between the point values is high {Cij ~ 1), as is the case for neighboring points of a fine lattice, one 
can develop series expansions of P n (u; Cij) in powers of (1 — Cij). We follow the calculations presented in Hamilton 
ct al. (1986), where a method for computing the series expansions for the distributions P2, P3, and P4 was presented, 
the expansion parameter being the lattice step A. As we have seen in Sec. Ill, knowledge of these three distrib ution s 
suffices to find the three Minkowski functionals on the plane. However, we do not assume an expansion of type ( 4.14 ) 
for the correlation function and keep (1 — Cij) as expansion parameters instead of the lattice step A, so that our results 
are applicable not only to smooth Gaussian fields. For completeness, we explain the method of calculation in some 
detail; this will also enable us to clarify the limitations of this method when applied to non-Gaussian distributions. 

We start with the Fourier representation of a Gaussian two-point symmetric probability density, 



P2 (x;C) 



d 2 k 

(2k) 2 



exp 



ikx- 



v' kx ~- 



P2 (k;C) 



Chk 2 



(Bl) 



The distribution (Bl) is symmetric with respect to an interchange of x%,2 and is characterized by the mean value 
(xi) = 0, variance {x 2 ) = 1, and correlation (x\X2) = C with \C\ < 1. The probability integral P 2 (u; C) is defined by 



P2 («; C) 



dx / dyp 2 (x;C) 



(B2) 
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This coincides with the distribution P 2 (it; C (A)) defined by Eq. (3.7). 
One then takes the partial derivative of P 2 (u; C) with respect to C . Since 

d . 



— p 2 (k;C) = -hk 2 p 2 (k;C), 



it follows that 



f) r°° r°° r ^v- 

— P 2 (u; C) = - J dxj dyj —^ kl k 2 e^p 2 (k; C) . 



(B3) 



(B4) 



To change the order of integration in Eq. (B4), it is necessary to replace the infinite upper limits in the integrations 
over x and y by a finite quantity M and take the limit of M — > 00 afterwards. The integrations over x and y yield 



M 



A I 



dx I dye lkx = 



iki 



ik-2 



The factor kik 2 in Eq. (B4) cancels. The distribution density p 2 (x; C) decays at spatial infinity, so 
-?-P 2 (u-C)= Urn \p 2 (u,u;C)-2p 2 (u,M;C)+p 2 (M,M ] C)}=p 2 (u,u;C). 



(B5) 



(B6) 



By means of this relation, the double probability integral in Eq. (B2) is reduced to a single integral 



P 2 (u; C) = P 2 («; 1) - / p 2 (tt, u; C) dC. 



(B7) 



Here we choose the boundary condition at C = 1 which corresponds to the limit of two coincident points, where the 
two-point distribution degenerates into a one-point Gaussian, 



P2 {x,y; 1) 



5{x-y) 



2tt 



ex P — h 



so it easily follows that 



F 2 (u; 1) = Pi («) = -( 1 - erf 



V2 



(B8) 



(B9) 



Now substituting into Eq. (B7) the explicit form of the distribution from Eq. (Al), 



p 2 (u,u; C) 



2ttVT - C 2 



exp 



1 + C 



(BIO) 



gives 



P 2 (u;C)=P 1 (u) 



dC 



c 



2nVl - C 2 



exp 



1 + C 



(Bll) 



Although the resulting integral cannot be evaluated analytically, it is easily expanded in powers of (1 — C) = e, 



f/.i: 



: exp 



1 + 3C 



2eexp 



1 - u 2 3 - 6u 2 + u 4 



12 



15 - 45u 2 + 15u 4 - m 6 3 105 - 420u 2 + 210w 4 - 28u 6 + u 
2688 6 + 55296 



160 

-e 4 + O (e 5 ) 



(B12) 



Another possible choice of the boundary condition is C = 0. The two-point distribution p 2 degenerates into a 
product of two one-point distributions pi, which leads to the following result, 



P 2 (u;C) = {P 1 (u)} 2 + [ 
Jo 



dC 
27iVl - C 2 



exp 



1 + C 



(B13) 
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As has been already noted by Hamilton et al. (1986), an extension of this technique to non-Gaussian distributions 
is problematic. Namely, one might consider the Fourier representation of a non-Gaussian two-point distribution 



p 2 (x; C, D, ...) = f d k 2 exp ikx-^ (kf + fcf) - Ck x k 2 - ^ / ), ,;/,■,/■•,/■•/ + 

J (27TJ L " 



(B14) 



where the coefficients Diji and the (omi tted ) higher-order terms of the series in k represent deviations from Gaussianity. 
Then a formal relation similar to Eq. (|B6|) may be derived: 



d_ 

dC 



P 2 (u;C,D,...)=p 2 (u,u;C,D,...). 



(B15) 



However, the limit of p 2 at C = 1 and fixed values of Diji and higher moments does not correspond to a limit of 
the two-point distribution with coincident points, since in that limit not only the second moment C, but also the 
higher moments D^i, ... are changed to make the expression under the exponential a function only of (ki + k 2 ) . 
Since p 2 (C = 1, D, ...) does not degenerate into a one-point distribution as before, one cannot make use of Eq. (B7) 
to compute the integral distribution P 2 . 

Having covered the derivation of the two-point distribution P 2 in detail, we now sketch the more complicated cases 
of the distributions P 3 and P4. The three-point distribution is 



P 3 (u;C) 



d.r 



dy / dzp 3 (x,y,z;C) 



(B16) 



for a symmetric three-point distribution density p% corresponding to the correlation matrix of Eq. (AE). The Fourier 
image of ps (x, y, z; C) is 



p 3 (fci,/c2,fc 3 ;C) =exp 



1 



{k\ +kl + kf)-C (fcife + k 2 k 3 + fc 3 fci) 



(B17) 



(B18) 



Differentiation of P3 (u; C) with respect to C yields a relation similar to Eq. (|B(i|). 

d f°° 
—P 3 (u;C)=3J dxp 3 {u,u,x;C) . 

At C = 1, the distribution P 3 (u; C) reduces to the same one-point distribution Pi (u) as before, so 

P 3 (u;C) = P 1 (u)-3 dC dxp 3 (u,u,x;C) . (B19) 

JO Ju 



The integral over x in Eq. (B19) can be evaluated by treating p 3 (u, u, x\ C) = f± (x) as a one-dimensional Gaussian 
distri bution in x (although /1 is not properly normalized). Its Fourier image /1 (fc) is easily found by integrating Eq. 
(B17) over k\ and k 2 : 



h (k) 



dk\dk 2 c ik lU +ik 2 u~ 
(2nf 



: exp 



p 3 {ki,k 2l k]C) 

1 r (1 - C) (2C + 1) k 2 2uC 



exp 



l + C 



2ttV1 - C 2 ~' r V l + C, 
For a general one-point Gaussian density given by its Fourier image 

pi (fc) = AT exp (~ak 2 /2 - i(3k) 

one readily obtains 



l + C 



dxpi (x) 



N 



1-erf 



u- j3 
V2a 



(B20) 



(B21) 



(B22) 



We only need to substitute the coefficients N, a and (3 from Eq. (B20) into Eqs. (B1S)-(B2C) to get an expression for 
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P 3 (u; C) = P 1 («) - 3 



dx 



: exp 




1 - erf ■ 



v /2(l + a;)(l + 2x) 



(B23) 



As before, the integrand can be expanded in powers of (1 — C) and integrated term by term. 

In a similar manner it is possible to evaluate the distribution P3 (u;Cij) for a general 3-point configuration with 
unequal correlations Cij between points. In that case, introduce a formal dependence of Cij on a parameter e, 



C,y(e) = l-e(l- dj) 



(B24) 



which interpolates between the original configuration (e = 1) and the degenerat e cas e Cy = 1 at e = 0. Then 
differentiate P 3 (u; Cij (e)) with respect to e to obtain a relation analogous to Eq. ( B18 ), and subsequently integrate 

over e from to 1. 

The calculation of the Euler characteristic on a hexagonal lattice as given by Eq. ( 3.39| ) can be further simplified 
because it contains only the combination Pi — 3P 2 + 2P 3 . With the distribution P 2 given by Eq. (B7), and using the 
identity 



Vi (u,u;C) 



dxp 3 (u, u, x; C) . 



one obtains 



Pi - 3P 2 + 2P 3 = 3 



dC 



dxp 3 (u, u, x; C) . 



In the notation of Eq. (B22) 

and we arrive at the expression 

Pi - 3P 2 + 2P 3 



-2 



dxpi (x) = iVerf 



u-/3 



2a 



dx 



27r\/r 



: exp 



erf ■ 



Uy/l 



+ xj yj2{l + x){l + 2x) 



(B25) 



(B26) 



(B27) 



(B28) 



After carrying out the series expansion in (1 — C) and term-by-term integration above, we obtained the result repre- 



sented in Eq. (4.21) of Sec. IV 



Finally we turn to the calculation of the distribution P4 (u). We only treat the special case when only two of the 
pair correlations between four points are independent. Namely, we assume that 



(0:10:3) = (X1X4) = (x 2 x 3 ) = (x 2 x 4 ) = B, 
(xix 2 ) = (x 3 x 4 ) = C. 

Then the Fourier image p 4 (k; B, C) of the 4-point distribution density p 4 (x; B, C) is 



p 4 (k; B,C) = exp 



-- (kf + kl + kl + kf) -B(ki + k 2 ) (k 3 + k 4 ) - C (kik 2 + k 3 k 4 ) 



(B29) 



(B30) 



The probability integral P 4 (it; B, C) is defined by Eq. (3J5). We shall obtain an expansion of P 4 in (1 — C) which will 
be useful for the case when the two pairs (0:1,0:2) and (0:3,0:4) of close-by points are separated by a comparatively 
large distance. Differentiating with respect to C gives the relation 



d_ 

dC J 



P 4 (u;B,C) — 2 / dx 4 / dx 2 p 4 (x 4 ,x 2 ,u,u; B,C) . 



(B31) 



As be fore, we choose the boundary condition at C — 1 where P 4 (u; B, 1) = P 2 (u ; B), and then we can integrate Eq. 
(B31) over C. The remaining task is to evaluate the double integral in Eq. (B31), and for this we again employ the 
same technique of converting two integrations into a single one. We treat the density p 4 [x 4 , x 2 , u, u; B,C) = f {x 4 , x 2 ) 
as a two-point (unnor maliz ed) Gaussian distribution for x\ and x%. The Fourier image of / {x\, x 2 ) is easily obtained 
by integration of Eq. ( B30| ), 
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/<'•■' • /'•->> = / [ ^e iuk *p i (k 1 ,k 2 ,k 3 ,k A ;B,C) 



N exp 



- o ( fc i + fc 2) - /Sfcifc - »7 (*i + *a) 



where 



iV : 



2vrVl -B 2 



exp 



1 + B 



a = 1 — 



2g 2 
1 + C" 



/3 = C- 



_2^2_ 

1 + C" 



2Bu 
1 + C 



(B32) 



(B33) 



variance 



The two-point density -/V" 1 / (xi, £2; ct, f3, 7) describes a Gaussian distribution with nonzero mean (xt)- 

(xf) = a, and correlation (2:1X2) — 0, and is similar to p 2 (x,y;C). Then we can directly use Eqs. ( Bll )-( Bl^ ) in 
which we must replace C by j3/a and u by (it — 7) /a/c* and also multiply by N: 



dx 



dy f (x, y; a, 0, 7) = — (l - erf ^J= 



Ndx 



: exp 



0-7) 2 

a(l+x) 



(B34) 



The resulting series in e = (1 — (3 /a) can then be converted into a series in (1 — C) , substituted into Eq. (B31) and 
integrated term by term over C . Because the resulting expressions arc quite cumbersome, we do not write them out 
explicitly. 

Similarly, one can obtain an expansion of P4 (u) for the case when all four points are close to each other and 
the parameters B and C approach 1 simultaneously; such an expansion is necessary for the calculation of the Eulcr 
characteristic on a square lattice. By parametrizing B = 1 — 5 and C = 1 — q8, where q is a formal constant parameter 
with value of order 1, and differentiating P4 with respect to d, one obtains a relation similar to Eq. (B31), 



d6 



Pa 0; B, C) 



dxi 



dx 2 [2qp4 (xi,x 2 ,u, u; B, C) + 4p 4 (xi,u,x 2 ,u;B, C)] . 



(B35) 



The double integrals are treated in the same manner as above, expanded in S and then integrated term by term. The 
actual values of B and C can be then substituted into the resulting series in S. In this case, however, the series (B12) 
cannot be used directly because in the limit of S = the parameters a, /3, 7 all tend to zero while the expansion 
parameter e = (1 — (3/a) of Eq. (B12) does not beco me small. Instead, the parameters a, f3, 7 can be expressed in 
terms of 6 and then a direct expansion of Eq. (B34) in S can be performed; the integrations for each term of that 
expansion can be calculated analytically. Again, we do write out the explicit form of the distribution P4. 
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